Bayesian segmentation of hyperspectral images 
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Abstract. In this paper we consider the problem of joint segmentation of hyperspectral images in 
the Bayesian framework. The proposed approach is based on a Hidden Markov Modeling (HMM) 
of the images with common segmentation, or equivalently with common hidden classification label 
variables which is modeled by a Potts Markov Random Field. We introduce an appropriate Markov 
Chain Monte Carlo (MCMC) algorithm to implement the method and show some simulation results. 

INTRODUCTION 

The most significant recent advances in remote sensing has been the development of 
hyperspectral sensors and software to analyze the resulting image data. Over the past 
decade hyperspectral image analysis has matured into one of the most powerful and 
fastest growing technologies in the field of remote sensing. The "hyper" in hyperspectral 
means "over" as in "too many" and refers to the large number of measured wavelength 
bands. Hyperspectral images are spectrally overdetermined, which means that they pro- 
vide ample spectral information to identify and distinguish spectrally unique materials. 
Hyperspectral imagery provides the potential for more accurate and detailed information 
extraction than possible with any other type of remotely sensed data. However the huge 
amount of data in hyperspectral images make its information exploitation difficult and 
image processing tools (classification, segmentation, comprising and coding) are needed 
to summarize the information included in these data. 

This paper will introduce a segmentation method for hyperspectral images. Several 
unsupervised and supervised algorithms have been developed for segmentation of mul- 
tispectral images. However, these algorithms fail to deliver high accuracies for classify- 
ing hyperspectral images Q3 El K |5]|. In this paper we consider the problem of joint 
segmentation of hyperspectral images in the Bayesian framework. 

The proposed approach is based on a Hidden Markov Modeling (HMM) of the 
images with common segmentation, or equivalently with common hidden classification 
label variables which is modeled by a Potts Markov Random Field. We introduce an 
appropriate Markov Chain Monte Carlo (MCMC) algorithm to implement the method 
and show some simulation results. This approach has previously been considered by [0 
for multi spectral images. 

In that work, the pixels of the same region in different images are assumed inde- 
pendent. This independence assumption is a valid hypothesis for multispectral images. 



However in hyperspectral images the pixel values in each channel are not independent. 
This work is then an extension to that work by considering a Markov model for these 
pixels along each channel. 

This paper is organized as follows: In the next two sections first we introduce our 
method for segmentation of hyperspectral images in the Bayesian framework. Then we 
propose an appropriate MCMC Gibbs sampling particularly designed for this segmen- 
tation task. Finally, in the last section we present some simulation results to show the 
performances of the proposed methods. 



MODELING FOR BAYESIAN SEGMENTATION 

Let gi(r) be the observed value of the pixel r, r £ Z 2 , in the spectral band i of a 
hyperspectral image. We model the observations by 

9i(r) = Mr)+ei(r), i = l,...,M (1) 

where /j(r) is the unknown perfect value of gi(r) and £j(r) is a noise. Note that if we 
consider images, the pixels r belong to a finite lattice S, and we will note S the number 
of pixels of this lattice. In the following we also use the notations 

9i = fi + £i or £=/ + e, (2) 

where g i = {gi(r),r £ S] and g = {g^i = 1,...,M} and a similar definition for / 
and e. We introduce a label variable z(r) for the regions and consider the region labels 
as common feature between all images. Thus the hidden variable z = {z(r),r £ S} 
represent a common classification of the images for different bands. The main result of 
this paper is estimation of joint segmentation label z. 

Assuming independent noises £j among the different observations we have 

M M 

p(g\f,e t ) = l[p(9i\fi) = UpeM-fi). (3) 

i=i i=i 

Assuming centered, white and Gaussian ~ A/"(0, cr| J), and S the number of pixels 
of an image, we have: 

p(^/ l ,^)=( 27r ^-) 2 exp|- 2 ^-||^-/ l || 2 }, (4) 

where cr| ~ XQ{oiq\Pq *) with unknown fixed parameters cuf i and [3^ (inverse gamma 
is conjugate prior for a random variance in the Gaussian case). 

To assign p(f i \z, .) we first define the sets of pixels which are in the same class: 

R k = {r:z(r) = k}, \R k \=n k 
f ih = {fr(r):z(r) = k}. 



We assume that all the pixels f ik of an image f i which are in the same class k will be 
Gaussian with a random mean m ik and a random variance af k , i.e. 

fi(r)\z(r) = k,m ik ,a 2 ik ~ Af(m ik ,a 2 k ) Vr G 5, (5) 



With these notations we have : 



P(fik\ m ik,^h) = [ / ~ j expl-^Vll/ife-^fc 1 !! 2 !- ( 6 ) 



and thus for i = 1, . . . , M 



K 



p(fi\z,m ik ,af k ) = II [ / H ^pl-^Vll/ife-^fciirK 

fe =i \y/2irat k J { Z(T ik ) 

where 1 is a vector with all components equal to 1, a 2 k ~ TQ(oti Q ,fli ), with unknown 
fixed parameters, and m ik is an autoregressive of order 1, AR(1), for each class k i.e. 

mk = 4>kmi-i k + Vik, O) 
where rj ik ~ AA(0, cr 2 2 ), 4> k and cr 2 are unknown fixed parameters. Therefore 

mi k \mi_ lk ~ A/"(0 fc mi_ lfc ,cr 2 o ). 

The assumption of ([7]) is the main difference of this paper with [6], i.e. in hyperspectral 
images the pixel values of a class k, in each channel, are not independent. 

Using the relation ([!]) and the density p(f i \z J m ik ,af k ) and p{£j), we can calculate 

9i\z,m ih ,a^ k ,al ~ M{m ik) a 2 ik + (T 2 ). (8) 

Finally we have to assign P(z). As we introduced the hidden variable z for finding 
statistically homogeneous regions in images, it is natural to define a spatial dependency 
on these labels. The simplest model to account for this desired local spatial dependency 
is a Potts Markov Random Field model: 



P(z) = -^exp<ja£ £ 5(z(r)-z(s))\, (9) 

1 \ a ) I r&SseV(r) 



where S is the set of pixels, 5(0) = 1, 5(t) = if t ^ 0, V(r) denotes the neighborhood of 
the pixel r (here we consider a neighborhood of 4 pixels), T(a) is the partition function 
or the normalization constant and a represents the degree of the spatial dependency of 
the variable z. 



ESTIMATION USING MCMC 



Let rrii = (m ik )k=i, k an d °f — {p^k)k=i,...,K be the means and the variances of the 
pixels in different regions of the images / i . We define L as the set of all the parameters 
which must be estimated in the Bayesian framework: 

O i = {o 3 ei ,m i ,(r*), i = l,...,M (10) 

and we note = {0i)i=i,...,M- Now we can write the expression of the joint posterior 
of p(f,z,0\g) by using the relations in the previous section. Then we propose the 
following Gibbs sampler, 

/ ~ f\g,z,e 

z ~ z\0,g 
~ 0\f,g,z 

to generate samples (f,z,0)^\ (f,z,0)( 2 \ - • -, and use them to compute any statistics 
(such as mean or median). We may note that in each of the previous Gibbs sampling 
steps, we again use Gibbs scheme to sample. For example f\g, z,0 by alternate sampling 
of f i |gfj, z, 0i. This procedure is also valid for 0. 

It can be shown that fi{r)\gi(r),z(r),&i has a Gaussian distribution and it can be 
sampled very easily. On the other hand, 

P(0i\fi,9i,z) ocpO^I/iiflO P(rni,<ri\fi,z). (11) 

For the last term p(mi,crf\f i ,z) we have to use a Gibbs algorithm and then sam- 
ple following the conditional distributions p(mj|<r?, f i7 z) and p(af\mi, f i ,z). It 
can be shown that <J%.\f i1 g i and ofj/^z have inverse gamma distributions and 

mk\fi,z,Vik> m i-ik~N(n 








V ik 


( n k 



(12) 

Note that if m ik s are independent as it is the case of multispectral images in H then we 
had 

Vik = v- k Mr + -2- E f^ r] 

\ a i0 a ikreR k 

where m i0 is a fixed number. 

Finally, we can write the posterior probability of z by 

M 

P(z\g,0)^p(g\z,0) P(z) = Y[p(9i\z,0i) P(z), (13) 



where can be calculated by using ([8]) and Q. As we choose a Potts Markov Random 
Field model for the labels z, we may note that an exact sampling of the a posteriori 
distribution P(z\g,0) is impossible. But we can still use a Gibbs sampling to generate 
parallel samples of z. 

For simplicity sake, we estimate the parameters (fi k , k = 1,...,K with a classical 
method and we consider it as constant in this section. If the series {mi fe }j = i v .. 5 M has 
an AR model, (k is fixed), then we can estimate <\> k efficiently, because the number of 
images M is large. 

Here we give the summary of the proposed algorithm for estimating z which has the 
following steps: 

1 . Find an initial joint segmentation of the hyperspectral image by using any simple 
segmentation or classification method such as k-means method, 

2. Calculate {mik}i=i,...,M for k = 1, . . . , K, 

3. Fit an AR model for each series {m ik } i= i^^ M , k — l,...,K with a classical 
method, 

4. Use the proposed Gibbs algorithm to generate samples of {m ik } using ( 12) and z 
using ([T3~]). 



SIMULATION 
Synthetic data 

To measure the performance of the proposed method, first we generate artificially a 
set of data, starting by a known segmentation z and generate the images fAz = k as 
homogeneous Gaussian with known mean m ik and variance af k with m ik = <p k mi-i k + 
r] ik , where \<p k \ < 1 is a fixed value. Then we generate the g { by adding a Gaussian noise. 
Finally we use these data as input for different segmentation methods and compare their 
relative performance. 

In this example, we know the original z. Therefore, comparison can be done by 
counting the number of misclassified pixels. 

Figure 1 shows the results of this experiment. The first row of Figure 1 shows z, f, 
and g. The second row shows the results of segmentation of images g by two methods: 
the first assume independence of m ik (i = 1, • • • , M) and second which is the proposed 
method use an AR modeling of m ik . The number of miss classified pixels for Zi is more 
than 1000 and for z 2 is less than 200. 




Z\ with independence and z 2 , with AR assumption on {m; fc } 

FIGURE 1. Bayesian segmentation of a simulated hyperspectral image with independence and AR 
assumption on {m^} with 50 iterations. 



Real data 



In the next step we applied the proposed method to a part of real data (Aviris). In 
Figure 2 we illustrate a real example hyperspectral images g which are 56, (128 x 124) 
images. This figure shows the reconstruction and joint segmentation results i.e. Z\ with 
independence and z 2 , with AT? assumption on {rrii k } after 200 iterations. Unfortunately, 
in this case we cannot give a quantitative comparing of the results. 




g and their marginals 




Z\ with independence and z 2 , with AR assumption on {m^} 

FIGURE 2. Bayesian segmentation of a part of real hyperspectral image with dimension 128 x 128 x 56 
with independence and AR assumption on {m 1( .} with 200 iterations. 
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